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Abstract We study size effects in the fracture strength of notched disor- 
dered samples using numerical simulations of lattice models for fracture. In 
particular, we consider the random fuse model, the random spring model 
and the random beam model, which all give similar results. These allow 
us to establish and understand the crossover between a regime controlled 
by disorder-induced statistical effects and a stress-concentration controlled 
regime ruled by fracture mechanics. The crossover is described by a scaling 
law that accounts for the presence of fracture process zone which we quan- 
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tify by averaging over several disordered configurations of the model. The 
models allow to study the development of the fracture process zone as the 
load is increased and to express this in terms of crack resistance (R-curve). 



1 Introduction 

Understanding how materials break is a fundamental open problem of sci- 
ence and engineering. The difficulties stem from the non-trivial dependence 
of the fracture strength on the characteristic lengthscales of the samples, as 
already noted by Leonardo da Vinci, who measured the carrying-capacity 
of metal wires of varying length (1). He observed that the longer the wire, 
the less weight it could sustain. The reason for this behavior is rooted in 
the structural disorder present in the material: the strength is dominated 
by the weakest part (subvolume) of the sample and its distribution can 
in principle be computed using extreme value statistics Longer wires 
have more weak parts and are thus bound to fail at smaller loads on aver- 
age. The quantitative understanding of this statistical size effect is difficult, 
since the important low-strength tails of the strength probability distri- 
bution are not easy to sample and since the material properties are often 
history-dependent. For instance, in quasi-brittle materials such as concrete 
and many other composites, where sample failure is preceded by significant 
damage accumulation ([31). 

An important engineering scenario and a typical experimental setting 
is to study the size effect in the presence of a pre-existing flaw, a notch. 
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Failure in this case is determined by the competition between deterministic 
effects, due to the stress enhancement created by notch, and the response of 
the disordered material around the defect to the stress concentration (|4|;|5|). 
This includes the stochastic damage accumulation. The effect of disorder 
can be treated in an effective medium sense by defining a Fracture Process 
Zone (FPZ) around the crack tip. For quasi-brittle materials, the size of 
this FPZ may not be negligible compared to the system size. Conversely, 
for small notches failure is more influenced by the FPZ than the notch 
itself, and depends on statistical disorder effects. Experimentally, it has 
been demonstrated that the critical crack may nucleate and propagate far 
from the pre-existing notch (|6|) in that case. 



The existing theories on the size effect start from linear elastic fracture 
mechanics (LEFM). Several formulations have been proposed in the liter- 
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2J). The 



ature and partly compared with experiments (|5|; [jj |8|; |8 ; 
problem is how to extend LEFM in the presence of disorder and concomitant 
damage. 



In LEFM the Griffith's stability criterion or equation for the maximum 
stress the specimen can bear reads a c ~ K c /^/aQ. Here ao is the linear 



size of the crack and the critical stress intensity factor K c ~ y/EG c is a 
function of the fracture toughness G c and the elastic modulus E (11 31). I n 



the size-effect law proposed by Bazant for quasi-brittle materials (5; 9; 



10|) 



the Griffith expression is generalized considering an additional lengthscale 
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£ due to the presence of a FPZ 



o c = K c /y/£ + a Q . (1) 

Equation ([T]) takes into account the both limits of a very large notch com- 
pared to the FPZ size and that of a very small notch. In the limit £/ao <C 
one has an expression that follows the LEFM scaling, a c ~ 1/^/ao- In the 
opposite limit of ao — > 0, the average strength is taken to be constant. 
Eq. |T|) has been shown to be in agreement with several experimental data 



sets (jlOj). However, three fundamental questions can be asked: first, does 
Eq. (fTJ) incorporate all the important effects? Second, what is the fracture 
toughness G c in the presence of disorder? Third, where does the FPZ scale 
£ originate and how does it depend again on the disorder? 

In this work we clarify the role of the disorder in the failure of notched 
quasi-brittle specimens using extensive simulations of disordered lattice 



models for fracture (14). A brief account of the results for the RFM has 
been published in Ref. (|15l ). In more detail, we consider the random fuse 
model (RFM), the random spring model (RSM) and the random beam 
model (RBM). We consider the failure of notched disordered samples and 
provide a microscopic justification of Eq. (p}. Studying the size scaling of 
strength by extensive numerical simulations is a difficult task due to the 
different length scales involved and the need of significant statistical aver- 
aging. We vary the disorder, which we model as a locally varying random 
failure threshold, and show that it plays a crucial role in determining the 
size effect. In particular, a lengthscale £ emerges from the simulations and 
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can be shown to be directly related to the FPZ size. Finally, for notch sizes 
smaller than a critical length a c , we observe a cross-over to the inherent, 
sample-size dependent strength of the unnotchcd sample. We outline the 
scaling theory to account for these effects. We show that the RFM results 
are confirmed in the RSM and the RBM. We also study the growth of the 
FPZ, showing that it is independent of the notch size ao and the system 
size L. 

2 Models 

2.1 Random fuse model 



In the RFM (|l6l), we consider a triangular lattice of linear size L with 
a central notch of length ao- Each fuse has the same conductance and a 
random breaking threshold t. This represents a locally varying fracture 
toughness/strength. The t lie between and 1, with a cumulative distri- 
bution P{t) = t 1 / , where D represents a quantitative measure of disorder. 
The larger D is, the stronger the disorder. The burning of a fuse occurs 
irreversibly, whenever the electrical current in the fuse exceeds breaking 
threshold t of the fuse. Periodic boundary conditions are imposed in the 
horizontal direction to simulate an infinite system and a constant voltage 
difference, V, is applied between the top and the bottom of lattice system 
bus bars. Numerically, a unit voltage difference, V = 1, is set between the 
bus bars and the Kirchhoff equations are solved to determine the current 
flowing in each of the fuses. Subsequently, for each fuse j, the ratio between 
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the current ij and the breaking threshold tj is evaluated, and the bond j c 
having the largest value, max., I 2 -, is irreversibly removed (burnt). The cur- 
rent is redistributed instantaneously after a fuse is burnt implying that the 
current relaxation in the lattice system is much faster than the breaking of 
a fuse. Each time a fuse is burnt, it is necessary to re-calculate the current 
redistribution in the lattice to determine the subsequent breaking of a bond. 
The process of breaking of a bond, one at a time, is repeated until the lattice 
system fails completely. In the present simulations, we have considered var- 
ious notch sizes for L = 64, 128, 192, 256, 320 and D = 0.1, 0.3, 0.5, 0.6, 0.75. 



2.2 Random spring model 



In the RSM, we consider a triangular lattice with nodes connected by linear 



springs with unit stiffness (jjJZ 
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2l|). As for the RFM, the bond 



breaking thresholds, t, are randomly distributed based on a thresholds cu- 
mulative probability distribution, P(t) = t 1 / for t G [0, 1]. The bond breaks 
irreversibly, whenever the force in the spring exceeds the breaking threshold 
force value, t, of the spring. Periodic boundary conditions are imposed in 
the horizontal direction and a displacement difference is applied between 
the top and the bottom of the lattice system. Numerically, a unit displace- 
ment, A — 1, is applied at the top of the lattice system and the equilibrium 
equations are solved to determine the force in each of the springs. Subse- 
quently, for each bond j, the ratio between the force fj and the breaking 
threshold tj is evaluated, and the bond j c having the largest value, maxj^-, 
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is irreversibly removed. The forces are redistributed instantaneously after a 
bond is broken implying that the stress relaxation in the lattice system is 
much faster than the breaking of a bond. Each time a bond is broken, it 
is necessary to re-equilibrate the lattice system in order to determine the 
subsequent breaking of a bond. The process of breaking of a bond, one at 
a time, is repeated until the lattice system falls apart. For the RSM, we 
consider a triangular lattice network of size L with a notch of size oo- In the 
present simulations, we have considered various notch sizes for L = 128, 256 
and D = 0.5,0.6. 



2.3 Random beam model 



231 ). we consider a two- 



In the random thresholds beam model (RBM) ([22]; 
dimensional triangle lattice system of linear size L. The vectorial RBM has 
three degrees of freedom (x-translation u, y-translation v, and a rotation 
8 about z axis) at each of the lattice nodes (sites), and each of the bonds 
(beams) in the lattice connects two nearest neighbor nodes. We assume that 
the beams are connected rigidly at each of the nodes such that the angle 
between any two beams connected at a node remains unaltered during the 
deformation process. These nodal displacements and rotations introduce 
conjugate forces and bending moments in the beam members. Using Tim- 



oshenko beam theory (|24r) . which includes shear deformations of the beam 
cross-section in addition to the usual axial deformation of cross-sections, 
the local stiffness matrix for a beam element that relates the local nodal 
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displacements and rotations to local nodal forces and bending moments in 
the beam's local coordinate system is given by 



K 
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where -E is the Young's modulus, G is the shear modulus, A is the beam 
cross-sectional area, / is the moment of inertia of beam cross-section, if, 
is the length of the beam, and a = ^rfp is the shear correction fac- 
tor, which denotes the ratio of bending stiffness to the shear stiffness. If 
shear deformation of beam cross-section is negligible, then a = and the 
Timoshenko beam theory reduces to Euler-Bernoulli beam theory. Equa- 
tion [2] presents a relation between local nodal displacements and rota- 
tions di = (uii,vu,8ii,uij,vij,8ij) T and local forces and moments F^ = 
(Fu,Vii, Mu, Fij ,Vij , Mij) T . In this setting, the subscript I refers to local 
quantities, the superscript T represents transpose of a vector or a matrix, i 
and j refer to i-th and j-th nodes of the beam, and F, V , and M refer to 
axial force, shear force, and bending moments respectively. 

Equilibration of the lattice system is achieved by first transforming these 
local quantities (d^ and F^) into global quantities d = (ui, Vi, 8i, Uj, Vj, 8j) T 



and F = (Fi,Vi, Mi, Fj,Vj, Mj) T through a coordinate transformation T 
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such that dg — Td, = TF, and K = T T K (oca ;T, and then satisfying 
equilibrium equations at each node such that 



£<ij>F x — 0; S < ij > F y — 0; S < ij > M — 



(3) 



where E<ij> implies that the summation is carried over all the intact bonds 
< ij > joining at node i. In the above discussion, the transformation matrix 
T is given by 

Q 



T = 



Q 



(4) 



where 



c s 

Q = -s c (5) 
1 

and c = cos(/3), s = sin(/3) refer to the direction cosines of the beam with 
(3 representing the angle between the beam axis and the cc-direction. 

In the present simulation, we start with a notched lattice system with 
beams having unit length, unit square cross-section and Young's modulus 
E = 1. This results in a unit axial stiffness [EAjl^ = 1) and bending stiff- 
ness (12EI/t\ = 1) for each of the beams in the lattice system. Since the 
beam can deform in two independent deformation modes (axial and bend- 
ing), we assume randomly distributed bond breaking axial and bending 
thresholds, t a and %, based on thresholds cumulative probability distribu- 
tions, P a (t a ) and Pb(h) respectively. As in the other models, the cumulative 
distributions arc defined as P(t) = t x l D in [0, 1]. 
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The failure criterion for a beam is denned through an axial force and 
bending moment interaction equation (similar to von-Mises criterion in 
metal plasticity) given by 



t a J h 

The beam breaks irreversibly, whenever the failure criterion r > 1. Periodic 
boundary conditions are imposed in the horizontal direction and a constant 
unit displacement difference is applied between the top and the bottom of 
lattice system. 

Numerically, a unit displacement, A = 1, is applied at the top of the 
lattice system and the equilibrium equations (Eq. [3]) are solved to determine 
the force in each of the springs. Solution of Eq. [3] results in global displace- 
ments and rotations d, using which the local displacements = Td and 
the local forces = Ki oca ;d^ are computed for each of the intact beams. 
Subsequently, for each bond k with nodes i and j, the quantities a k = (j^j 
and bk — max (l^'l'l M J D are evaluated, and the bond k c having the smallest 
value, 



-b k + y/tij + 4a k 

is irreversibly removed (When = 0, then rk = ^-). The forces are re- 
distributed instantaneously after a bond is broken implying as in the other 
models that the stress relaxation in the lattice system is much faster than 
the breaking of a bond. Each time a bond is broken, it is necessary to re- 
equilibrate the lattice system in order to determine the subsequent breaking 
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Fig. 1 The stress concentration profiles for the RFM (left), the RSM (center) 
and the RBM (right). 

of a bond. The process of breaking of a bond, one at a time, is repeated 
until the lattice system falls apart. 

3 Strength and size effects 

We perform numerical simulations of the models discussed above, concen- 
trating on the failure strength. Notice that the three model differ mainly in 
the way stress is redistributed. To illustrate this point, we report in Fig.[T]the 
stress concentration profiles in a triangular lattice (L = 512) with a notch of 
size ao = 16. Although the angular distribution of the stress profiles differs, 
the way stress decay from the crack tip is very similar, approaching for large 
L the 1 / yfr decay expected from the theory of elasticity. 

In order to obtain reliable results, the strength should be averaged over 
several realizations of the disorder. In the present simulations, we have used 
a minimum of N r = 200 realizations and in some cases up to N r = 8000 
realizations. Fig. [5^ reports the strength cr c , averaged over different con- 
figurations, with varying ao, D, and L for the RFM. The most instructive 
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way of plotting is to consider the inverted square strength, X/cq. Assuming 
Eq. ([1]), \fa\ should become a linear function of ao for large enough notches. 
Plotting the data in this way in Fig.[2}D shows that for ao ^> 1, the scaling of 
Eq. (fTJ) is recovered asymptotically. Extrapolating the linear part towards 
ao = 0, we can define a disorder-dependent intercept £(-D), that should be 
related to the FPZ size. Furthermore, the slope of the linear part of the data 
(1/K^(D)) is also disorder-dependent, which implies a disorder-dependent 
fracture toughness G C (D). Finally, a careful observation reveals that for 
small ao less than a critical crack size a c , the strength scaling crosses over 
from a stress concentration dominated LEFM scaling [Eq. fl}] to a disorder 
dominated scaling (see Fig. [3]). That is, for ao <§C a c , the strength scaling 
deviates significantly from Eq. fl} and saturates to a value that depends 
on disorder and the sample size, <j c (L, D). In particular, the strength of the 
unnotched system (for ao = 0) is finite and is smaller than the LEFM limit 
ifc/i/f given by Eq. fl]). In Fig. [3J we present a comparison between REM, 
RSM and the RBM. The general features of the strength are the same in the 
three models, indicating that only the decay of stress concentration 

is relevant for the size effect, while the precise angular dependence of stress 
concentration around a notch is not important. 



In Ref. (|15l ) we presented a scaling theory that extends the earlier scal- 
ing law given by Eq. (TTJ) beyond its actual regime of validity. A correct 
scaling expression has to accommodate the three separate phenomena vis- 
ibile in Fig. [2] for small notches, the dominance of statistical effects that 
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dictate a c (L,D), the cross-over to the LEFM-like regime, and then finally 
an Equation ([1]) like scaling at large ao- 

The cross-over takes place at a scale a c above which a c follows Eq. ([1}. 
For small notches, ao <C a c the strength is determined by extremal statis- 



tics as is appropriate in the limit ao — > ()14l ). Then one expects to see a 
weak size effect, typically logarithmic in L. In real materials, the scaling will 
depend on the damage accumulation and the defect populations that exist 
in the specimens. cr c (L, D) is not a constant however, as Eq. ([1]) would pre- 
dict, and deviates significantly from the LEFM-based theory, which would 
in general predict that the samples arc weaker than their actual strength 
a c (L,D). 

The location of the cross-over (notch size) a c follows by equating the 
strength prediction of Eq. JT]) and the scaling of notchless specimes, l/a(L,D) 2 
(a c + S)/K"l. An appropriate scaling theory, valid for all a , is then 

= £ + a /(a c /a ) (8) 
where the statistical physics -like scaling function f(y) has the limits 

!1 if y < 1 
(9) 
y if v > l 

Thus we have for the cross-over scale 

a c ~(K c (D)/a c (L,D)) 2 -£(D). (10) 

For ao > a c , fracture is governed by LEFM and a scaling of the kind of 
Eq. ([T]) is recovered. The effect of disorder, according to the scaling theory, 
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is incorporated in the three parameters a c (L, D)), K C (D), and £,(D). In the 
following, we first discuss the first two parameters and return to £ below 
in much more detail. Qualitatively (since the behaviour of a c (L,D) is an 
independent issue entirely), the effect of changing disorder strength for a 
fixed L can be seen as follows. For stronger disorder, the cross-over will 
take place at larger a c , since the stress concentration of the notch will be 
screened (as we demonstrate below). At fixed disorder, a c increases with 
L since notchless specimens get weaker. The fracture toughness G c (since 
in the models E = 1) seems in our simulations to be proportional to the 
average model element strength at weak disorder at least, and perhaps gets 



reduced with strong disorder (|15f). More numerical work in this direction 
might be interesting. 



4 The fracture process zone 

Our numerical simulations allow to monitor the damage evolution prior to 
failure and can thus be used to study the development of the FPZ. For 
a single realization of the disorder, we only see diffuse damage up to the 
peak load, and it is difficult to determine the size of the FPZ. On the 
other hand, the FPZ can be clearly measured after averaging the damage 
over several realizations of the disorder. Hence, the size of FPZ should thus 
be considered in statistical terms as the region around the crack tip where 
damage is most likely to occur. Considering for simplification a projection of 
the average damage along the notch main axis direction, we obtain a profile 
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Fig. 2 a): The strength in the RFM for several disorders D, notch sizes ao and 
system sizes L. b): a scaling plot of the data according to Eq. |T}. 
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Fig. 3 A scaling plot of the strength according to Eq. |T} for the RFM (for 
L = 192, D = 0.6), the RSM (L = 256 D — 0.5) and the RBM (L = 256 D = 1.0) 
models. The qualitative features of the data for different system sizes and disorders 
are the same in all the three models. 



that is decaying exponentially towards a homogeneous background value: 
d(x) — A + B exp (— 2x/£fpz) (see Fig. 0]). The factor 2 in the exponential 
comes from the fact that in our geometry the FPZ extends from the two 
edges of the notch. We have analyzed the data for different values of do and 
L in order to check that the profiles do not depend on L and on an, as long 
as this is not too far from a c . For an << a c we naturally do not expect to see 
such a "damage cloud" around the original defect. However, it seems likely 
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that one could measure £ around the most critical microcrack. Recall that 
in this regime one expects the strength to saturate at a c = \J K c /(a c + £). 

Notice that the LEFM stress intensity factor would indicate a l/\/r - 
like divergence of the stress at the crack tip. It is evident that the observed 
exponential decay of damage profile d is in contrast to a 1/s/r -like decay 
and should be naturally interpreted as a screening of the crack tip caused by 
the disorder. In fact, the FPZ size £,fpz depends on the disorder strength 
D as shown in Fig. [5l The data can be roughly described by a power law 
relation £,fpz ~ D 3 ^ 2 . As discussed in Ref. (|15l ). if we plot the fracture 
process size £,fpz{D) against the intrinsic scale £ resulting from the fits of 
the strength data to Eq. (pj, we obtain a linear relation. Hence, we can 
conclude that £ is indeed a direct measure of the FPZ size. Fig. [5] reports a 
comparison of the damage profiles in RFM, RSM and RBM. It can be seen 
that the results are qualitatively similar for all the three models considered. 
The value of the FPZ size £fpz, however, differs slightly for the for the 
three models. 

The FPZ progressively develops before the peak load by damage accu- 
mulation. To visualize this process, we have computed damage profiles at 
different values of the applied stress. One can then obtain the FPZ size 
£,fpz{D) as a function of the stress. As can be seen in Fig. [7J there is a 
gradual increase of £,fpz{D) with stress. This growth relates to the R-curve 
of the material Jlli which is usually defined in terms of the elastic energy 
released due to crack growth G = dU /da (jj). In the RFM model, we can 
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derive G from the lattice "elastic" energy U = I 2 /(2S), where S is the 
conductivity, as 

AU _ I 2 AS 

where AS — (So — S) is the conductivity change after the crack has ex- 
tended by Aa such that a = ciq + Aa, and So is the initial conductivity. 
We report the R-curve for different values of disorder in Fig. [5] The data 
is shown for the current values I E [I c /2,] in which £(J) can be extracted 
from the damage profiles with a reliable accuracy. In this regime the R- 
curves show in general two behaviors: an initial rapid increase due to the 
accumulating damage that changes the average conductivity, followed by a 
saturation as the FPZ starts to increase even more rapidly. The conductivity 
change is a mean-field phenomenon that accounts for the total damage in 
the system, proportional to D(I), and thus to the failure thresholds P(i c ). 
The growth of £(/) is not expected to be so simply related to D. It is inter- 
esting to note that as a result the R-curves for various disorders overlap in 
the manner depicted in Fig. [8] The size-dependence of these R-curves would 
be expected to be negligible as long as the strength is governed by Eq. ([!} . 

5 Conclusion 

We have resorted to simulations of statistical fracture models to analyze the 
problem of the size-effect in the failure of quasi-brittle materials. For large 



notches, the simulations recover the expression based on LEFM Q_5|; |9 



10|) 



and allow to relate the effective FPZ size £ to the actual average damage 
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Fig. 4 a): Damage profiles along the crack axis for various notch sizes ao. Damage 
profiles follow an exponential decay on a uniform damage background, i.e., d(x) — 
A + Bexp (— 2x/£fpz), where A and B are constants and x is the distance from 
the crack tip along the crack axis. In order to show that £fpz is independent 
on ao, we report in the inset the profiles after subtracting the background and 
normalizing so that the curves superimpose, b) Damage profiles for two different 
lattice sizes L and two different notch sizes ao. The profiles do not depend on L. 




Fig. 5 The dependence of the FPZ size on disorder for the RFM (L = 128, 
ao = 16). Estimates from the strength (£) and from damage profiles (,fpz are 
similar. We could not obtain reliable estimates from damage profiles for very 
weak disorder (D = 0.1). 



profiles. As the notch size is decreased we observe a crossover at a novel 
scale a c to a disorder-dominated size-dependent regime that is not described 

n 

by LEFM and is furthermore seen in experiments (|15l). All the regimes are 
summarized in a generalized scaling expression for the strength of disordered 
media. 

Several interesting future questions remain, like theoretical computa- 
tions of parameters such as a c , K Cl and the detailed understanding of the 
origin and shape of the statistical FPZ. These would be in particular im- 
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Fig. 6 A comparison of the damage profiles measured in RFM, RSM and RBM 
using system size L = 128, disorder D — 0.6, and an initial notch size ao = 16. 
Two thousand samples are used for averaging the damage profiles. The result show 
that while the profiles are qualitatively similar ^fpz and the damage backgound 
differ. 



portant in order to help to achieve practical predictions. Recall our results 
have shown, that all such parameters are dependent on disorder, which in 
the models used translates into the damage accumulated at a given local 
stress. This quantity is evidently hard to access experimentally, but is un- 
fortunately theoretically necessary. For such reasons, it would be relevant 
to investigate three dimensional systems and possibly look at other kinds 
of disorder (eg. locally varying elastic moduli). 
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Fig. 7 The FPZ size as a function of the applied stress normalized by the peak 
stress as obtained from damage profiles. Data are for RFM for different values of 
disorder (oo =16, L = 128). 
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